Executive summary

Niniejszy raport zawiera analizę danych pacjentów chorych na koronawirusa przyjętych do szpitala Tongji w Wuhan na przełomie stycznia i lutego 2020 roku. W raporcie uwględnione zostały analizy wraz z wykresami przedstawiające zależności pomiędzy różnymi cechami pacjentów a śmiertelnością. Analizowanymi atrybutami były między innymi: dane demograficzne pacjentów (płeć, wiek, etc.) oraz dane medyczne pacjentów powstałe w wyniku wykonanych na nich badaniach krwi. W raporcie przedstawione zostały również korelacje pomiędzy różnymi biomarkerami oraz dołączone zostały wykresy zmienności najbardziej kluczowych atrybutów na przestrzeni kolejnych badań. W ostatnim rozdziale raportu stworzony został klasyfikator dokonujący predykcji śmiertelności pacjenta zakażonego koronawirusem. Dokładność opisanego klasyfikatora osiągnęła poziom 95.40%.

Wykorzystane biblioteki

library(readxl)
library(httr)
library(knitr)
library(plyr)
library(gridExtra)
library(tidyverse)
library(janitor)
library(scales)
library(kableExtra)
library(qwraps2)
library(rmarkdown)
library(caret)
library(reshape2)
library(plotly)
library(naniar)
library(gganimate)

Wczytywanie danych

Surowe dane zostały pobrane z adresu url (http://www.cs.put.poznan.pl/dbrzezinski/teaching/zed/wuhan_blood_sample_data_Jan_Feb_2020.xlsx) bezpośrednio bez konieczności ręcznego pobierania zbioru danych w celu zapewnienia powtarzalności eksperymentów.

GET(url, write_disk(data_file <- tempfile(fileext = ".xlsx")))
raw_df <- read_excel(data_file)

Przetwarzanie zbioru danych

Rozmiar danych

Zbiór danych składa się z 6120 rekordów medycznych opisywanych przez 81 atrybutów, z których pierwsze 7 odnosi się bezpośrednio do samego pacjenta (jego identyfikator, data przyjęcia do szpitala, wiek, …), a pozostałe 74 atrybuty zawierają informacje o wynikach przeprowadzonych badań.

Pacjenci

Zbiór danych zawiera informacje o 375 pacjentach, u których stwierdzono obecność choroby COVID-19. Dane zostały zebrane na przestrzeni od 2020-01-10 do 2020-03-04.

Atrybuty główne - informacje o pacjentach

  • patient_id (PATIENT_ID) - unikalny identyfikator przypisany do każdego pacjenta,
  • test_date (RE_DATE) - pole, którego znaczenie w dokumentacji zbioru nie jest jednoznacznie opisane, ale na podstawie pozostałych atrybutów zbioru można wyprowadzić, że jest to data wykonania badania
  • age - wiek pacjenta
  • gender - płeć pacjenta, przyjmuje wartości 1 lub 2, gdzie 1 oznacza mężczyznę, a 2 kobietę. Przypisane wartości do płci nie było jawnie udokumentowane, natomiast poprawną identyfikację umożliwiały dopiero opublikowane statystyki zbioru.
  • admission_date (Admission time) - data przyjęcia pacjenta do szpitala.
  • discharge_date (Discharge time) - data wypisania pacjenta ze szpitala w przypadku ozdrowienia lub data śmierci.
  • survived (outcome) - wynik przebytej przez pacjenta choroby, przyjmuje wartości 0 lub 1, gdzie 0 oznacza wyleczenie, a 1 - śmierć. Przypisanie odpowiedniego wyniku do wartości atrybutu było możlwie dzięki dostępnym statystykom zbioru, podobnie jak w przypadku atrybutu ‘gender’.

Dane demograficzne pacjentów

Wśród pacjentów najliczniejszą grupę pod względem płci stanowią mezczyzni - 224 pacjentów (59.7%).

Średnia wieku pacjentów wynosi 58.83, pełne informacje o rozkładzie wieku pacjetnów dostępne są w tabeli poniżej.
Min Median Mean Max
18 62 58.83 \(\pm\) 16.46 95

Dodatkowo do danych o pacjentach dodany został nowy atrybut wyliczeniowy - age_group przypisujący pacjentów do dziesięcioletnich grup wiekowych, powstały na podstawie podziału atrybutu numerycznego - age. Najliczniejszą grupą wiekową wśród pacjentów jest (60,70] - 114 pacjentów (30.4%), natomiast najmniej liczną grupą jest (10,20] - 2 pacjentów (0.5%).

Statystyki dotyczące demografii pacjentów zostały również zwizualizowane na wykresach kołowych poniżej:

Hospitalizacja

Minimalnym czasem hospitalizacji był 0:02:01:57 (dd:hh:mm:ss), a maksymalnym - 35:04:05:54. Średnim czasem jaki pacjenci spędzili w szpitalu był 10:20:29:08.

Pełna dystrybucja czasu hospitalizacji jest przedstawiona w tabeli poniżej:

age_group Mean Max Min Standard_deviation
(10,20] 8:16:05:49 14:15:13:49 2:16:57:48 8:10:25:16
(20,30] 11:13:04:00 17:16:22:40 0:10:44:51 4:22:31:21
(30,40] 14:02:16:38 35:04:05:54 0:03:50:46 7:02:59:42
(40,50] 13:13:26:15 29:04:45:21 0:12:57:09 7:12:09:37
(50,60] 11:07:50:50 33:17:10:13 0:03:40:29 7:16:55:57
(60,70] 10:10:54:04 29:00:38:10 0:03:55:53 7:13:14:12
(70,80] 7:18:15:24 27:18:53:29 0:11:46:57 6:14:10:13
(80,90] 9:00:41:43 31:16:54:28 0:02:01:57 8:04:20:22
(90,100] 4:10:55:24 4:18:31:13 3:21:22:44 0:09:38:49

Liczba przeprowadzonych badań

Pacjenci byli poddawani badaniom ze średnią częstotliwością 2.26 badania na dzień. Największą ilością przeprowadzonych badań podczas pobytu pacjenta w szpitalu było - 59, natomiast najmniejszą - 1.

Śmiertelność

Bezwględna śmiertelność pacjentów ze zbioru danych wyniosła 46.0%.

Tabela z procentowo wyrażoną śmiertelnością w zależności od grupy wiekowej
survived (10,20] (20,30] (30,40] (40,50] (50,60] (60,70] (70,80] (80,90] (90,100] Combined
yes 50.00% 100.00% 97.73% 81.08% 65.15% 44.44% 12.96% 14.81% 0.00% 54.02%
no 50.00% 0.00% 2.27% 18.92% 34.85% 55.56% 87.04% 85.19% 100.00% 45.98%

Najbardziej narażeni byli pacjenci z grupy wiekowej (90,100], których śmiertelność wyniosła 100.0%, ale jest to również jedna z najmniej licznych grup pacjentów.

Najmniej narażoną grupą byli pacjenci w wieku (20,30], u których nie zanotowano przypadków śmiertelnych.

Dodatkowo ze względu na małą reprezentację grupy wiekowej (10,20], tylko dwóch pacjentów - przypadek śmiertelny oraz ozdrowiciel, współczynnik śmiertelności (oraz przeżywalności) jest przekłamany.

Biomarkery

Pozostałe atrybuty zbioru danych (znajdujące się w kolumnach od 8 do 81) zawierają wyniki badań pacjentów. Podstawowe statystyki dotyczące biomarkerów znajdują się w tabeli poniżej:

Częstotliwość mierzonych czynników medycznych podczas badań pacjentów

Najczęściej mierzonymi czynnikami podczas przeprowadzonych badań były: liczba leukocytów (White blood cell count) oraz liczba erytrocytów (Red blood cell count).

Najrzadziej sprawdzane czynniki medyczne wsród pacjentów
Czynnik Liczba badan % badan
Interleukin 10 267 4.36%
Interleukin 2 receptor 268 4.38%
Interleukin 8 268 4.38%
Tumor necrosis factora 268 4.38%
Interleukin 1ß 268 4.38%
Interleukin 6 272 4.44%
Najczesciej sprawdzane czynniki medyczne wsród pacjentów
Czynnik Liczba badan % badan
Red blood cell count 1127 18.42%
White blood cell count 1127 18.42%
Serum potassium 980 16.01%
calcium 979 16.00%
serum sodium 975 15.93%
Serum chloride 975 15.93%

Korelacje atrybutów

Poniższa macierz zawiera wartości korelacji pomiędzy atrybutami biomedycznymi testowanych wśród pacjentów. Istnieje kilka korelacji pomiędzy atrybutami których wartość wynosi 1. Dzieje się tak ze względu na dwa czynniki: duża liczbę przeprowadzonych badań i mała oraz zmienna liczba sprawdzanych czynników podczas pojedynczego badania. Przykładowo korelacja czynników hemoglobin oraz HCV antibody quantification wynosi 1, ale liczba rekordów w jakich te dwa czynniki były wspólnie mierzone wynosi 2, przez co wartość ich korelacji może być błędnie zawyżona.

Najważniejsze atrybuty biomedyczne

W artykule powiązanym z analizowanym zbiorem danych - Yan, L., Zhang, HT., Goncalves, J. et al. An interpretable mortality prediction model for COVID-19 patients - przedstawionych zostało 10 najważniejszych atrybutów zbioru danych na podstawie wyników modelu klasyfikatora Multi-tree XGBoost. Trzema najważniejszymi wyselekcjonowanymi cechami były (w kolejności malejącej ważności atrybutów): Lactate dehydrogenase (dehydrogenaza mleczanowa, LDH), lymphocyte count (liczba limfocytów) i High sensitivity C-reactive protein (wysoko czułe białko C-reaktywne, hs-CRP).

Poniżej przedstawiony został wykres, pochodzący z wymienionego artykułu, przedstawiający 10 najbardziej kluczowych atrybutów zbioru danych przy szacowaniu śmierelności pacjentów.

Autorzy m.in. poniższych publikacji:

potwierdzają obserwacje, iż badanie poziomów: enzymu dehydrogenazy mleczanowej (LDH), wysokoczułego białka C-reaktywnego (hs-CRP) oraz liczby limfocytów pacjenta mogą służyć przy szacowaniu śmierelności/dotkliwości pojedynczych przypadków koronawirusa.

Dehydrogenaza mleczanowa (LDH)

Wysokoczułe białko C-reaktywne (hs-CRP)

Limfocyty

Klasyfikacja

Na podstawie opisanej wyżej oceny ważności atrybutów, wybranymi atrybutami do przygotowania modelu klasyfikacji zostały: Lactate dehydrogenase (dehydrogenaza mleczanowa, LDH), lymphocyte count (liczba limfocytów) i High sensitivity C-reactive protein (wysoko czułe białko C-reaktywne, hs-CRP).

Czyszczenie i przygotowanie zbioru danych do klasyfikacji

W zbiorze danych występują uszkodzone rekordy - są to wiersze, w których nie ma żadnych wyników badań pacjentów oraz brakuje daty wykonania badania (kolumna test_date, RE_DATE w nieprzetworzonym zbiorze danych). Rekordy uszkodzone zostały usunięte w przetwożonym zbiorze danych.

Liczba rekordów przed usunięciem - 6120, liczba rekordów po usunięciu - 6106.

Ze względu na dużą liczbę badań oraz na to, że nie każde badanie zawierało taki sam zestaw mierzonych czynników zbiór danych został pogrupowany ze względu na pacjentów, a wartości parametrów: LDH, limfocyty i hs-CRP zostały wypełnione średnimi zmierzonymi wartościami niniejszych atrybutów dla poszczególnych pacjentów.

Podział zbioru danych

Zbiór danych został podzielony w proporcji 75% danych przydzielonych do zbioru uczącego (parametr p), natomiast podanym atrybutem do stratyfikacji (parametr y) był atrybut decyzyjny survived.

inTraining <- createDataPartition(
        y = patients_results$survived,
        p = .75,
        list = FALSE)

Model klasyfikacyjny

Dane zostały poddane powtarzanej kroswalidacji oraz podany został zakres parametru mtry którego wartości będą służyły do optymalizowania klasyfikatora.

rfGrid <- expand.grid(mtry = 10:30)

gridCtrl <- trainControl(
    method = "repeatedcv",
    summaryFunction = twoClassSummary,
    classProbs = TRUE,
    number = 10,
    repeats = 5)

Wybranym modelem klasyfikatora wykorzystywanym do predykcji śmiertelności pacjentów został Random Forest. Wartość parametru mtry=29 została dobrana automatycznie podczas procesu uczenia.

set.seed(23)
fitTune <- train(survived ~ .,
             data = training,
             method = "rf",
             metric = "ROC",
             preProc = c("center", "scale"),
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 10)

fitTune
## Random Forest 
## 
## 264 samples
##   4 predictor
##   2 classes: 'no', 'yes' 
## 
## Pre-processing: centered (4), scaled (4) 
## Resampling: Cross-Validated (10 fold, repeated 5 times) 
## Summary of sample sizes: 237, 237, 238, 238, 238, 238, ... 
## Resampling results across tuning parameters:
## 
##   mtry  ROC        Sens       Spec     
##   10    0.9759524  0.9366667  0.9456190
##   11    0.9744603  0.9333333  0.9512381
##   12    0.9735992  0.9333333  0.9500952
##   13    0.9753849  0.9383333  0.9484762
##   14    0.9757579  0.9333333  0.9458095
##   15    0.9756548  0.9350000  0.9540000
##   16    0.9745833  0.9400000  0.9459048
##   17    0.9747381  0.9400000  0.9457143
##   18    0.9745397  0.9400000  0.9470476
##   19    0.9760992  0.9400000  0.9454286
##   20    0.9762698  0.9383333  0.9470476
##   21    0.9762976  0.9316667  0.9499048
##   22    0.9731310  0.9266667  0.9441905
##   23    0.9742063  0.9366667  0.9526667
##   24    0.9730635  0.9316667  0.9482857
##   25    0.9746667  0.9333333  0.9498095
##   26    0.9773135  0.9366667  0.9526667
##   27    0.9727619  0.9383333  0.9495238
##   28    0.9765516  0.9433333  0.9513333
##   29    0.9766071  0.9383333  0.9512381
##   30    0.9741667  0.9316667  0.9484762
## 
## ROC was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 26.

Predykcja

Macierz pomyłek

Poniżej przedstawiona została macierz pomyłek wygenerowana dla wyuczonego klasyfikatora. Miara false positive dla tego klasyfikatora wynosi - 0.05, a miara false negative wynosi - 0.9787234

no yes
no 38 2
yes 1 46

Miary klasyfikatora

Dokładność klasyfikatora wynosi - 96.55%, a natomiast miara Kappa przyjmuje wartość - 93.05%.

Szczegółowe wartości przedstawione są w tabeli poniżej:
score
Sensitivity 0.9743590
Specificity 0.9583333
Pos Pred Value 0.9500000
Neg Pred Value 0.9787234
Precision 0.9500000
Recall 0.9743590
F1 0.9620253
Prevalence 0.4482759
Detection Rate 0.4367816
Detection Prevalence 0.4597701
Balanced Accuracy 0.9663462